#include "Quad8NodeShpFunc.h"

using SG::DataStructure::Common::Real;
using SG::Algebra::Matrixd;
using SG::FEMSolver::Element::Quad8NodeShpFunc;

Quad8NodeShpFunc::Quad8NodeShpFunc ()
:NSet (), DNDxiSet ()
{}

/** @brief 计算积分点处形函数及形函数对参数坐标导数
 * 
 * 缓存各个积分点处的计算结果，避免下一次积分点循环时重复计算（内存换取计算时间）
 */

void Quad8NodeShpFunc::Compute (const SG::DataStructure::FEM::IntegPointArray& points)
{
    const auto nGauss = points.size ();
    NSet.resize (nGauss);
    DNDxiSet.resize (nGauss);
    const Real quarter{0.25};
    const Real half{0.5};
    for (std::size_t index = 0; index < nGauss; ++index)
    {
        const Real& xi{points[index].m_x};
        const Real& eta{points[index].m_y};

        NSet[index] = Matrixd (8, 1);
        NSet[index] (4, 0) = half * (1.0 - xi*xi) * (1.0 - eta);
        NSet[index] (5, 0) = half * (1.0 + xi) * (1.0 - eta * eta);
        NSet[index] (6, 0) = half * (1.0 - xi*xi) * (1.0 + eta);
        NSet[index] (7, 0) = half * (1.0 - xi) * (1.0 - eta * eta);
        NSet[index] (0, 0) = quarter * (1.0 - xi) * (1.0 - eta) - half * (NSet[index] (4, 0) + NSet[index] (7, 0));
        NSet[index] (1, 0) = quarter * (1.0 + xi) * (1.0 - eta) - half * (NSet[index] (4, 0) + NSet[index] (5, 0));
        NSet[index] (2, 0) = quarter * (1.0 + xi) * (1.0 + eta) - half * (NSet[index] (5, 0) + NSet[index] (6, 0));
        NSet[index] (3, 0) = quarter * (1.0 - xi) * (1.0 + eta) - half * (NSet[index] (6, 0) + NSet[index] (7, 0));

        auto DNDxi = Matrixd (2, 8);
        //边中节点处
        DNDxi (0, 4) =xi * (eta - 1.0);
        DNDxi (1, 4) =half * (xi * xi - 1.0);

        DNDxi (0, 5) = half * (1.0 - eta * eta);
        DNDxi (1, 5) = -1.0 * eta * (1.0 + xi);

        DNDxi (0, 6) = -1.0 * xi * (1.0 + eta);
        DNDxi (1, 6) = half * (1.0 - xi * xi);

        DNDxi (0, 7) = half * (eta * eta - 1.0);
        DNDxi (1, 7) = eta * (xi - 1.0);

        // 四个顶点
        DNDxi (0, 0) = quarter * (-1.0 + eta) - half * (DNDxi (0, 4) + DNDxi (0, 7));
        DNDxi (1, 0) = quarter * (-1.0 + xi)  - half * (DNDxi (1, 4) + DNDxi (1, 7));
        DNDxi (0, 1) = quarter * ( 1.0 - eta) - half * (DNDxi (0, 4) + DNDxi (0, 5));
        DNDxi (1, 1) = quarter * (-1.0 - xi)  - half * (DNDxi (1, 4) + DNDxi (1, 5));
        DNDxi (0, 2) = quarter * ( 1.0 + eta) - half * (DNDxi (0, 5) + DNDxi (0, 6));
        DNDxi (1, 2) = quarter * ( 1.0 + xi)  - half * (DNDxi (1, 5) + DNDxi (1, 6));
        DNDxi (0, 3) = quarter * (-1.0 - eta) - half * (DNDxi (0, 6) + DNDxi (0, 7));
        DNDxi (1, 3) = quarter * ( 1.0 - xi)  - half * (DNDxi (1, 6) + DNDxi (1, 7));
        // 储存计算结果
        DNDxiSet[index] = DNDxi;
    }
}